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Abstract 

The  overall  goal  of  this  project  is  to  develop  a  new  and  efficient  boundary  element 
method  (BEM)  program  for  the  stress  analysis  of  3-D  anisotropic  rock  masses  that  are  perturbed 
by  irregular  topographies  and  underground  openings.  This  stress  analysis  will  provide  direct 
information  on  the  best  (or  optimal)  selection  of  the  surface  location  and  hit  angle  for  a 
penetrator.  We  have  identified  three  tasks  to  achieve  this  goal,  which  include  theoretical  and 
analytical  development;  computational  and  numerical  development;  and  laboratory  investigation 
and  field  validation.  This  final  technical  report  presents  our  accomplishments  related  to  the  first 
and  second  tasks,  that  is,  the  theoretical  and  analytical  development,  and  computational  and 
numerical  development. 

First,  the  analytical  solution  and  niimerical  implementation  are  presented  for  elastostatic 
displacement  Green’s  fimction  for  3D  rock  masses  of  general  anisotropy.  Excerpts  from  the 
authors’  FORTRAN  code  are  included.  A  numerical  algorithm  for  Ae  calculation  of  the 
derivatives  of  the  Green’s  displacements  and  stresses  is  also  introduced.  Secondly,  these  Green’s 
functions  are  incorporated  into  a  BEM  code  developed  by  the  authors.  Thirdly,  numerical  results 
of  Green’s  displacements,  stresses  and  stress  derivatives  are  presented  and  compared  to  the 
closed-form  solutions  for  transversely  isotropic  rocks.  Finally,  the  BEM  code  based  on  the 
current  Green’s  functions  is  tested  and  the  numerical  results  are  compared  to  those  using  a  BEM 
code  based  on  the  exact  Green’s  functions.  It  is  shown  that  the  Green’s  functions  derived  in  this 
report  are  accurate  and  the  corresponding  BEM  code  is  correct.  This  BEM  code  is  now  ready  for 
the  laboratory  comparison  and  field  validation. 
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1.  Introduction 

Green’s  functions  are  important  in  the  formulation  of  boundary  integral  equations  and  in 
the  solution  of  those  equations  by  the  boundary  element  method  (BEM)  [1,2].  In  order  to  handle 
rock  and  rock  mass  anisotropy,  the  static  Green’s  functions  in  a  3-D  anisotropic  full-space  are 
required.  Previously,  several  methods  were  proposed  to  calculate  these  Green’s  functions.  These 
include  the  numerical  integral  method,  series  expansion  technique,  dual  reciprocity  technique, 
and  the  eigenvalue/eigenfunction  method.  While  the  first  three  methods  are  approximate,  the  last 
one  requires  solving  a  6x6  eigenvalue  system. 

After  reviewing  thoroughly  this  topic,  Wang  [3]  derived  explicit  expressions  for  3D 
elastostatic  Green’s  displacements  in  general  anisotropic  solids  and  integrals  of  Green’s 
displacement  derivatives  over  segments  and  rectangles. 

At  the  outset,  this  report  reviews  some  of  the  basic  concepts  inherent  in  Wang’s 
formulation  for  Green’s  displacements.  A  particularized  account  of  the  authors’  derivation  and 
implementation  of  these  expressions  follows,  along  with  the  key  parts  of  the  authors’  own  code 
written  in  FORTRAN. 

Secondly,  a  numerical  algorithm  for  the  calculation  of  the  derivatives  of  the  Green’s 
displacements  and  stresses  is  subsequently  introduced;  it  allows  the  discretization  of  the 
boundary  to  be  of  the  most  general  type  in  a  BEM  formulation. 

Thirdly,  it  is  described  how  the  preceding  implementations  were  utilized  within  a  BEM 

code. 

Finally,  numerical  examples  of  Green’s  displacements,  stresses  and  stress  derivatives  are 
presented  for  a  transversely  isotropic  solid,  so  as  to  allow  a  comparison  with  a  previously 
available  closed-form  solution  [4].  Numerical  examples  of  BEM  calculations  are  also  given  and 
the  results  are  compared  with  exact  solutions  and  previously  published  munerical  results. 
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2.  Outline  of  the  analytical  solution 


2.1  Notation 


Consider  the  geometry  of  Figure  la  where  (O,  xi,  X2,  X3)  is  a  Cartesian  coordinate  system 
in  a  3-dimensional  Euclidean  space  R^,  (ui,  U2,  U3)  the  corresponding  right-handed  orthonormal 
basis  and  x=(xi,  X2^,jC3)  a  point  in  this  space.  We  ^sume  that  the  anisotropic  body  is  embedded  in 
this  space. 

Let  n=(ui,  ^2,  «3)  be  a  vector  whose  components  in  are  rii,  n2,  «3,  with  respect  to  uj,  U2, 
U3,  respectively.  We  can  imagine  n=(«i,  «2,  also  as  a  point  in  a  Cartesian  coordinate  system  of 
a  3-dimensional  Euclidean  space,  which  will  be  called  the  n-space.  Let  Q  be  any  closed  surface 
containing  the  origin  of  the  n-space  and  dQ(n)  an  infinitesimal  area  element  of  this  surface 
around  point  n=(ni,  W2,  ^3),  see  Figure  1  .b. 

Throughout  this  paper,  indicates  the  dot  product  of  two  vectors  and  “x”  indicates  the 
cross  product.  Also,  a  comma  indicates  partial  differentiation  with  respect  to  a  variable,  i.e. 

/i=—  and  summation  over  repeated  indices  is  assumed. 

dXi 


2.2  Basic  equations 

Consider  an  unbounded  homogeneous  anisotropic  linearly  elastic  solid  subjected  to  a 
point  load  in  the  fixed  coordinate  system  (O,  xi,  X2,  X3)  depicted  in  Figure  l.a.  The  Green’s 
function  will  be  denoted  by  gpk(x)  and  gives  the  displacement  in  the  Xp-direction  at  x  produced  by 
a  point  load  applied  at  the  origin  O  in  the  Xk-direction.  Let  <7,-,  be  the  stress  tensor,  Ui  the 
displacement  field  and 


2  ^ 


(la) 


the  infinitesimal  strain. 

The  constitutive  relation  of  linear  elasticity  reads: 

~  ^iJpy^pg 

where  cypq  is  the  elastic  tensor,  which  is  fully  symmetric  and  positive  definite,  i.e.: 

^ypg  ~  ^jipg  ~^tjqp  ~  ^pgy 


(lb) 

(lc) 


Cij-p^aijOpg  >0  V  non-zero  tensor  ay  (Id) 

Inserting  the  kinematics  relation  (la)  into  the  constitutive  relation  (lb)  and  the  latter  into 
the  equilibrium  equation: 

+  =0  (2a) 

where  Fi  is  the  body  force  per  unit  volume,  we  obtain  the  following  three  second-order  partial 
differential  equations: 
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(2b) 


C..  U 


+  Fi=0 


once  the  symmetry  (Ic)  of  the  elastic  tensor  is  taken  into  account.  Because  the  Green’s  function 
is  relative  to  a  point  force,  (2b)  becomes: 

^ijpqS pkjq^^  ~  (^) 


where  6ik  is  the  Kronecker  delta  and  6(x)  the  Dirac  delta. 

In  the  subsequent  implementation,  the  6x6  matrix  D  of  Voigt  constants  is  introduced  such 

that; 

* 

(Jk  =  D  8k  (4) 

where  a  kK^^nfo  <^22k,  crsak,  023k,  <yi2k)  is  the  Green’s  stress  vector  (relative  to  a  point  force 
applied  in  direction  k  at  the  origin)  and  Ek=(siik,  e22k,  Sssk,  2s23k,  28i3k,  2Ei2k)^  is  the  Green’s 
strain  vector  (relative  to  a  point  force  applied  in  direction  k  at  the  origin)  with  its  components 
being  defined  as : 


i^ik,J  S jk,i ) 


The  follovwng  rules  apply  between  the  components  Cyks  and  Dap  and  account  for  the  first 
two  symmetries  of  formula  (Ic)  [5]: 


If  i=J,  a=i ; 

If  a=9-i-j ; 


If  A=5,  ^=k ; 

If  k^s,  P=9-^-j' ; 


(6) 

(7) 

(8) 
(9) 


The  third  symmetry  of  formula  (Ic)  is  accomplished  by  the  symmetry  of  D,  therefore: 


If(X<p,  Cyks  Dap;  (10) 

If  a>p,  Cijks  =  Dpa;  (11) 


23  Analytic  solution 

We  notice  that  CijpqWjWq  is  symmetric  and  positive  definite,  so  that  it’s  inverse  is  well 
defined.  We  set: 

^ipM=^ijpqf^j^q  ’  r^W  =  (wA)’‘  (12) 

Consider  now  the  following  identities,  in  which  integration  is  taken  in  the  n-space  over 
any  closed  surface  Q  including  the  origin  (see  Figure  lb)  and  use  is  made  of  Eqs.  (A.6b)  and 
(A.  8)  (see  Appendix  A): 
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l5=n-x 


^  Jr;;J(nWn  •  X>/Q(n)  = 

'■  •^  Q  n 


=  (Cprgs^r^s  )"*  =  S^j, 

n  *="•*  Q 


f/Q(ii)  = 
c/Q(n)  = 


(13) 


n  PI 


Since  the  last  member  in  (13)  can  be  wrijyten  in  terms  of  the  plane  representation  of  the 
delta  function  (A.  10),  we  have: 


-  Jr;Kii)^(n  •  xVQ(n)  =  -U^Si^Six) 
n 


(14) 


By  comparing  (3)  with  (14),  we  get  the  following  integral  expression  for  the  Green’s 
displacement  function: 


Spk(^)  =  Jr;l (ii)5(n  •  x)dQ{n) 

n 


(15) 


We  now  try  to  write  (15)  in  a  form  suitable  for  integration  by  means  of  the  Theory  of 
Residues  [6].  To  this  end,  we  express  the  inverse  tensor  r;^(n)  as: 


r;i(n)= 


(16) 


where  Ap^{n)  is  the  adjoint  matrix  of  cypqWjnq  and  D(ii)  is  the  determinant  of  CijpqWjnq,  e.g.: 

^pM  =  a(^[cijpgnjng\-,  D{vL)  =  AQ\[cyp^njn^\  (17) 

It  is  convenient  to  use  coordinates  (^,  ti)  introduced  in  Appendix  B.  Using  Eq.  (B.4), 
equation  (15)  becomes: 


(18) 


where  r  is  the  distance  between  the  field  and  source  points  (see  equation  (B.l.a)). 

Since  rpjt(n)  is  a  3x3  matrix,  each  entry  of  its  adjoint  matrix  Apk  is  a  polynomial  of  order 

4  in  1^  and  q.  Similarly,  determinant  D(n)  is  a  polynomial  of  order  6  in  4.  C  and  q. 

Following  Wang  [3],  we  choose  £2  as  a  rectangular  parallelepiped  having  size  2Lx2Lx2 
(Fig.  2)  and  let  the  dimension  L  go  to  infinity.  Since  over  surfaces  other  than  Si  and  S2  the 
integrand  in  (18)  approaches  zero  as  Ml},  the  contribution  of  the  integration  over  every  surface 
but  Si  and  S2  vanishes.  Moreover,  the  integrand  in  (18)  is  symmetric  with  respect  to  since  only 
even  powers  of  ^  are  involved.  This  leads  to  twice  the  integration  over  Si  (^=1): 
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(19) 


1 

AtP' 

We  recall  that  [6],  for  any  infinitely  differentiable  function  /x)  and  any  positive  real 
number  a: 


^d{ax)f{:)Cpx 

—00 


—  [8(x)f(xla)dx  =  —f{(i) 
a  i  a 


(20) 


which  is  justified  by  the  change  of  variables  ax-»x.  Therefore,  integrating  (19)  with  respect  to  r\ 
yields : 


(21) 


The  integral  in  (21)  can  be  handled  by  means  of  the  Theory  of  the  Residues  [6].  The  poles 
are  the  roots  of  the  polynomial  equation  of  sixth  order  in  (^: 

Dip +  ^q)  =  0  (22) 


If  11560  is  real,  rpjt(n)  is  a  real  symmetric,  positive  definite  matrix  and  then  its  eigenvalues 

are  all  real  and  positive  [6],  thus  its  determinant  can  never  be  zero.  If  we  want  Eq.  (22)  to  be 
satisfied,  ^  must  therefore  be  complex.  A  corollary  of  the  Fundamental  Theorem  of  Algebra  [7] 
tells  us  that  a  real  polynomial  of  order  N  has  N  roots  and  that  if  Cf=a+ib  is  a  root,  its  conjugate 
C*=a-ib  must  also  be  a  root  [7].  In  our  case  we  have  3  roots  satisfying  : 

f>(p+4'mq)=o  (23) 


with: 

lm^„>0  /m=1,2,3  (24) 

and  we  may  write; 

iXp  +  fq)  =  ‘  -  f.  Xf  -  ♦)  (25) 

it=0  m=l 

where  ntk  are  the  coefficients  of  the  sextic  polynomial  £)(p  +  ^q)  with  respect  to  C,.  Equation  (23) 
is  called  the  sextic  equation  of  elasticity  and  it  has  been  shown  that  no  closed-form  solution 
exists  for  its  roots  [8].  Thus  this  equation  must  be  solved  munerically.  Integral  (21)  can  now  be 
expressed  in  terms  of  the  residues  at  the  poles,  taking  into  account  that  it  must  be  real; 


4*(p-^Cq) 


07(C  -i.  -Ck  •) 


k=\ 

k^m 


(27) 
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It  is  the  authors’  experience  that  the  6  poles  resulted  simple  even  with  isotropic  materials. 
Should  the  poles  be  multiple,  a  slight  change  in  the  elastic  constants  will  result  in  single  poles, 
with  negligible  errors  in  the  computed  Green’s  tensor. 

Some  features  of  formula  (27)  need  to  be  pointed  out: 

1)  since  rp;t(n)  is  symmetric,  its  adjoint  matrix  Apk  is  also  symmetric  and  so  is  Green’s  tensor 

gpk.  As  a  consequence,  only  6  terms  out  of  9  must  be  calculated; 

2)  for  two  points  xi  and  X2  aligned  along  the  same  line  passing  through  the  origin,  the  summation 

over  index  m  has  the  same  value; 

3)  as  a  consequence  of  2),  gpk  approaches  zero  as  1/r  when  r->oo; 

4)  as  a  consequeiii5e  of  2),  g-pk  depends  only  on  ttie  relative  position  of  the  source  point  and  field 

point.  Thus  the  implementation  can  proceed  considering  the  source  point  always  at  the  origin, 
by  an  applicable  translation;  this  leads  to  an  important  simplification  of  the  implementation 
itself; 

5)  the  numerical  solution  of  a  polynomial  of  sixth  order  is  the  only  numerical  step  required  in 

order  to  obtain  the  entire  Green’s  function. 

3.  Implementation  of  the  analytic  solution 

The  key  steps  in  the  implementation  of  the  analytic  formulation  fall  essentially  into  two 

groups: 

1)  entries  of  (p  +  ^q)  in  terms  of  stiffness  matrix  components  Dap,  coordinates  of  field  point  x, 

and  variable 

2)  coefficients  Oi  i=l,...,7  of  sextic  polynomial  (25). 

In  order  to  keep  the  expressions  as  simple  as  possible,  vector  v  introduced  in  Appendix  B 
must  coincide  with  one  of  the  base  vectors  ui,  ua,  U3.  In  the  following,  we  have  assumed  that  v  is 
either  (1,0,0)  or  (0,1,0);  the  choice  between  them  must  be  done  on  the  basis  of  x  and  affects  step 
1)  only. 

3.1  Entries  of  matrix  (p  +  ^q) 

Each  entry  of  matrix  (p  +  ^q)  is  a  polynomial  of  order  two  in  ^  of  the  form: 

Tj/  (p  +  Cq)  =  byi  +  by2C  +  (28) 

As  an  example,  coefficients  iyk  are  given  in  Appendix  C  in  FORTRAN  format  for  the 
case  v=(l,0,0);  they  take  into  accoxmt  Eqs.  (6)-(ll)  so  that  they  are  functions  of  the  entries  of 
matrix  D,  which  are  much  more  manageable  then  the  fourth-order  tensor  cypq.  Coefficients  feyk 
are  also  functions  of  the  field  point  coordinates.  The  correspondence  between  notation  used  in 
the  text  and  notation  used  in  the  Fortran  code  is  established  in  Table  1.  Mathematica  2.2  was 
used  to  get  the  expression  of  coefficients  ^yk. 

3.2  Coefficients  of  the  sextic  polynomial 

Although  it  is  possible  to  get  the  expressions  of  coefficients  Ui  (/=!,. ..,7)  directly  in  terms 
of  the  stiffiiess  matrix  components  Dap  and  of  the  field  point  coordinates,  this  leads  to 
expressions  as  long  as  40  pages,  which  cannot  be  handled  easily  by  the  compiler  and  are  not 
computationally  efficient.  Thus  the  idea  is  to  derive  them  in  terms  of  the  6yk  coefficients 
introduced  in  (28).  Determinant  D(p  +  i^q)  can  in  fact  be  written  as  the  sum  of  trinomials: 


7 


"^13^22^31 +2ri2r23r3i  -riir23r32  ~ri2r2ir33 +rj|r22r33 


(29) 


If  we  put  (28)  into  (29),  we  realize  that  each  coefficient  a\  (relative  to  C*)  is  a  sum  of  trinomials 
of  6ijk  coefficients.  Let  /,  m,  n  be  the  last  indexes  of  the  coefficients  of  a  trinomial,  we  found  that 
these  indexes  have  to  satisfy  I  +  m  +  n=i-2  in  order  for  the  product  to  be  of  order  i-1.  For 
example,  only  coefficients  biju  such  that  k=l  contribute  to  ai.  Thus,  in  general : 


/,7n,«:/+m+n=/— 2  /,m,w:/+/n+«=/-2  /,m,f*:/+m+«=/-2 


^12f^2 1 jhsJc  +  ^  ^  1(^22  j^33i  4, 

/,w,«:/+m+w=/-2  /,m,/i:/+m+n=j‘“2 


(30) 


Coefficients  atj  calculated  according  to  (30)  are  given  in  Appendix  D  in  Fortran  format 
(see  Table  1  for  notation  correspondence).  It  is  to  be  noted  that  the  present  implementation  of  the 
analytic  solution  is  quite  simple,  when  compared  to  the  complexity  of  the  problem  in  hand. 


4.  Derivatives  of  the  Green’s  displacements  and  stresses 

Analytic  solutions  for  the  integral  over  segments  and  rectangles  of  the  Green’s 
displacement  derivatives  were  proposed  by  Wang  [3].  However,  expressions  of  the  derivatives  of 
the  Green’s  displacements  are  necessary  if  the  boundary  discretization  must  be  more  general. 
Several  attempts  were  made  by  the  authors  in  order  to  get  a  closed  form  solution  of  the  Green’s 
function  derivatives.  The  most  promising  of  them  started  from  the  derivation  of  (21)  and  led  to 
the  expression : 


/w=l 


-J-ImV 

2;z7'  ■“ 


‘-lie. -c.  ’Of.  -ft  ‘) 

k^m 


dx. 


k=l 

k^m 


+^lmy  Ito 

Inr  dx, 

m=\  ^  L  s. 


(31) 


In  this  equation,  the  last  term  is  very  complicated,  therefore,  a  numerical  algorithm  was 
implemented  based  on  the  Lagrange  polynomials  [9].  Despite  of  its  simplicity,  this  approach  has 
been  proven  to  be  efficient,  accurate,  and  robust. 

Following  [9],  let  a  function  J(x)  be  known  at  the  n  points  xi<  X2<...<  Xn.  Let  us  call 
yi=^Xi)  and  set : 


n  n  n 

A:=l  r=l  r=l 

r^k  r^k 


(32) 
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The  complete  Lagrange  interpolation  function  is  : 

n\ 


where : 


(33) 


(34) 


By  taking  the  deriy^tive  of  (33)  and  evaluating  it  in  Xr,  we  get : 

n\ 

Thus,  the  error  in  the  first  derivative  is  ; 

n\ 


(35) 


(36) 


F'(x^)  is  nothing  but  the  product  of  the  distances  between  Xr  and  the  other  chosen  abscissas. 
Therefore,  if  the  intervals  between  the  chosen  abscissas  is  constant,  its  minimum  value  is  attained 
at  the  mid  point  (or  two  mid-points  if  n  is  even)  of  the  segment  between  x\  and  Xn.  It  follows  that 
the  best  approximation  of  the  value  of  /'(xj  obtainable  using  the  polynomial  derivative  is 
attained  at  the  mid  point  (or  two  mid  points  if  n  is  even)  of  the  segment  between  x\  and  Xn. 

It  can  be  shown  that  the  derivative  of  the  Lagrange  polynomial  is  [9] : 


^  -  Xi 


Jt=l 

k^r 


Pr-yic 


fM 

P'M. 


(37) 


If  we  choose  a  polynomial  of  order  2,  i.e.  3  abscissas,  from  (35)  and  (37)  we  get: 
/'(^2)  =  ^[-/(^i)+/fe)]-|- (38) 


where  h  is  the  distance  between  two  consecutive  abscissas  and  ^2  is  a  point  comprised  between  xi 
and  X3. 

Now,  let  X  be  the  field  point  at  which  we  want  to  calculate  the  Green’s  stress  component 
Gijk  defined  in  Section  2.2.  To  this  end,  the  expression  of  the  Green’s  strain  component  Sijk  at  x  is 
necessary  in  order  to  calculate  aijk  by  means  of  (4).  Using  (38),  the  following  expressions  of  the 
Green’s  strain  component  syk  at  x  are  obtained: 


M..gu(x  +  Ai)-git(x-Ai) 

^  2/2 

(x) «  g2it(x  +  A2)-g2;^(x-A2) 


1h 


(39) 


(40) 


(41) 
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(42) 


^  M  +  A2)-  ^3Jt(x  -  ^2)+  glki'^  +  A3)-  glki^  -  A3) 
SmW - Yh 

^  L.\  ■  gl/t(x  +  A3  )-  gi;t(x  -  A3  )+  +  Al  )-  giki^  -  A,  ) 

ml  )~  2/2 

^  M  . .  gl;t(x  +  A;)-  gi^(x  -  A2)+  g2kix  +  Al)-  g2t(x  -  AJ 

^Vlk\  /  7 


(43) 

(44) 


where  Ai=(A,0,0),  A2=(0,  h,  0),  A3=(0,0,  h). 

In  order  to'^get  the  complete  Green’s  stress  and  strain  at  point  x,  it  is  thus  necessary  to 
compute  the  Green’s  tensor  at  6  points  in  the  neighborhood  of  x.  The  choice  of  interval  A  is  a 
crucial  decision.  An  extensive  numerical  investigation  has  led  us  to  the  conclusion  that  the  best 
value  of  the  interval  is: 


h  =  rl0-^ 


(45) 


where  r  is  the  distance  between  the  field  and  source  points. 

It  is  also  noteworthy  that  the  attempts  aimed  at  increasing  the  accuracy  of  the 
approximation  by  adding  other  terms  to  (39)-(44)  led  to  no  appreciable  improvement.  In  the 
authors’  opinion,  the  reason  for  the  good  performance  of  this  scheme  lies  primarily  in  the  smooth 
and  monotonic  behavior  of  the  Green’s  displacements. 

In  order  to  calculate  internal  stresses,  the  derivatives  of  the  Green’s  stresses  and 
displacements  are  needed  with  respect  to  the  coordinates  of  the  source  point  (see  Section  5).  By 
virtue  of  Observation  4  in  Section  2,  these  derivatives  are  equal  but  opposite  in  sign  to  the 
derivatives  taken  with  respect  to  the  field  point  Xi.  According  to  (38),  the  derivatives  of  the 
Green’s  stresses  are  approximated  as: 


2h 


/=1,2,3 


(46a) 


and  the  derivatives  of  the  Green’s  displacements  as: 


gg(x  +  A2)-g,;^.(x-A2) 

2h 


7=1, 2,3 


(46b) 


5.  BEM  formulation 

Consider  an  elastic  body  (finite  or  infinite)  with  the  following  displacement  and  traction 
conditions  imposed  on  the  boundary  r=ru+r: 

My(x)=M^.(x)  xeTu  (47a) 

xeTt  (47b) 

with  ttj  being  the  external  normal  to  Ft. 

For  each  internal  point  Xp  the  following  integral  equation  holds  [10]: 
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(48) 


r  r 

where  C/^(x^,x)  and  7],  (x^,x)  are  the  Green’s  displacements  and  tractions,  respectively.  By 
virtue  of  Observation  4  in  Section  2,  f7*-(x^,x)  and  7)y  (x^,x)  are  equal  to: 

(49) 

,  (50) 

m 

If  Xp  approaches  a  point  Xb  on  the  boundary,  (48)  is  modified  as: 

4"y(**)+  J^y(x6,x)w^(x>r(x)=  Jc/‘.(xj,x)7;.(x>r(x)  (51) 

r  r 

where  dy  are  coefficients  that  depend  only  on  the  local  geometry  of  the  boundary  at  Xb. 

The  term  on  the  right-hand  side  of  (51)  has  weak  singularity  (see  Observation  3,  Section 
2)  and  can  thus  be  integrated  by  means  of  a  usual  Gauss  quadrature  technique.  The  rigid-body 
motion  method  [11]  can  be  used  to  overcome  the  Cauchy-type  singularity  in  the  first  integrand 
and  at  the  same  time  to  avoid  the  calculation  of  coefficients  diy 

Eqimtion  (51)  can  be  discretized  and,  once  boimdary  conditions  (47)  are  taken  into 
accoimt,  the  resulting  algebraic  system  of  equations  can  be  solved  for  the  unloiown  boundary 
displacements  and  tractions.  Then  Eq.  (48)  can  be  used  to  calculate  internal  displacements. 

In  order  to  get  the  internal  stresses,  it  is  necessary  to  take  the  derivative  of  (48)  with 
respect  to  the  internal  coordinates  Xp.  This  yields: 


j7]*,/(x^,x)M/x>r(x)=  Ji7;_/(xp,x)r/x>r(x) 

r  r 

(52) 

where: 

1 

b 

1 

(53a) 

(53b) 

with  o-^.  ,(x  -  xj  given  by  (46a)  and  giji{x  - x^)  by  (46b). 

Once  «,,/(x^)  are  obtained,  the  internal  stresses  are  calculated  by  means 

of  the  following 

equation  similar  to  (4)  and  (5): 

ct  =  D  s 

(54a) 

where: 

<^“(<^11,  022,  cr33,  <^23,  ^13,  <^12)^ 

(54b) 

E  (Sll,  C22>  S33j  2E23,  26i3,  28i2) 

(54c) 

^iJ  ~  ^j,i ) 

(54d) 

11 


The  implementation  described  in  Sections  3  and  4  was  incorporated  into  an  existing  3D 
BEM  code  (see  [1 1]  for  its  description)  according  to  the  procedure  presented  in  this  Section. 

6.  Numerical  examples 

6.1  Green’s  displacements,  stresses  and  derivatives  of  the  stresses 

Let  us  consider  a  transversely  isotropic  and  linearly  elastic  solid,  whose  plane  of 
transverse  isotropy  is  parallel  to  the  xiX2  plane.  A  closed-form  solution  exists  in  this  case  [4]  for 
the  Green’s  displacements,  stresses  and  derivatives  of  the  stresses.  This  solution  will  be  used  (as 
implemented  in  [1 1])  to  validate  the  proposed  formulation. 

The  material  properties  are  as  follows:  .£=20-10^  kNW,  .£’=4-10'^  kN/m^,  v=0.25, 
v’=0.25,  G -1.6-10'’  kN/m^  where  E  and  are  the  Young’s  modulus  in  the  plane  of  transverse 
isotropy  and  in  the  direction  normal  to  it,  respectively;  v  and  v’  are  the  Poisson’s  ratios 
characterizing  the  lateral  strain  response  in  the  plane  of  transverse  isotropy  to  a  stress  acting 
parallel  and  normal  to  it,  respectively;  G  ’  is  the  shear  modulus  in  planes  normal  to  the  plane  of 
symmetry.  The  corresponding  stif&iess  matrix  D  is  (only  the  upper  half  is  given): 

^88  72  40  0 

88  40  0 

24  0 

16 

sym. 

\ 

The  field  point  is  placed  at  x=(- 1,0.8, 1.5)  m.  The  displacements,  stresses,  stress 
derivatives  are  collected  in  Tables  2,  3  and  4,  respectively.  The  agreement  between  the  closed- 
form  solution  and  the  present  formulation  is  very  good  for  all  three  quantities.  Only  for  some 
components  of  the  derivatives  of  the  stresses  (Table  4)  is  the  relative  difference  not  negligible; 
however,  it  must  be  noted  that  the  magnitude  of  these  components  is  small  with  respect  to  the 
remaining  components,  thus  leading  to  negligible  errors  in  the  computation  of  the  internal 
stresses,  as  will  be  shown  in  Section  6.2. 

6.2  BEM  models 

The  geometry  of  the  following  two  examples  is  the  same.  A  cube  of  transversely  isotropic 
material,  whose  edge  is  Im  long,  is  discretized  very  coarsely  with  6  nine-node  isoparametric 
elements  and  with  a  total  of  26  nodes  (see  Figure  3).  The  faces  of  the  block  are  parallel  to  the 
coordinate  planes.  The  block  is  subjected  to  uniform  compression  of  1  kNW  on  the  two  faces 
parallel  to  die  xiX2  plane.  For  this  case,  an  exact  solution  exists  [12]. 

In  the  first  example  we  adopt  the  same  material  constants  as  considered  in  Section  6.1. 
However,  the  plane  of  symmetry  is  no  longer  parallel  to  the  X\X2  plane,  but  is  inclined  with  dip 
orientation  cp=60°  and  dip  angle  \)/  =45°  (see  Figure  4).  Consequently,  stiffness  matrix  D  referred 
to  x\,  X2,  X3  axes  is  now  fully  populated  (only  the  upper  half  is  given): 


0  0^ 
0  0 
0  0 
0  0 
16  0 
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■10''kN/m^ 


(55) 
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(56) 


(14.5 


9.9 

9.6 

-0.2 

18.5 

9.6 

-2.2 

12.8 

-1.6 

3.2 

sym. 


-3.1177  -1.5588^ 
-1.0392  -1.9052 
-2.7713  0 

0  -1.0392 

3.2  -0.2 


•lO^kN/m^ 


This  case  was  also  considered  in  [1 1]  where  it  was  solved  using  the  closed-form  solution  for  the 
Green’s  displacenjgnts,  stresses  and  stress  deriv^ives  by  Pan  and  Chou  [4].  The  results  obtained 
with  these  formulations  and  with  the  present  one  are  given  in  Tables  5  and  6;  they  indicate  that 
for  both  the  displacements  (on  the  boimdary  and  internal)  and  internal  stresses,  the  exact  values 
and  the  values  calculated  in  [1 1]  are  in  very  good  agreement. 

The  second  example  resembles  the  one  reported  in  [13].  Here  a  (transversely  isotropic) 
zinc  cube  is  considered,  whose  plane  of  symmetry  is  parallel  to  the  X\X2  plane.  The  stiffness 
matrix  referred  to  xi,  X2,  x^  axes  is  now  (only  the  upper  half  is  given): 


^161.00  34.20  50.10  0 

161.00  50.10  0 

61.00  0 
38.3 


V 


0 

0 

0 

0 

38.30 


0 

0 

0 

0 

0 

63.40 


GPa 


(57) 


It  is  to  be  noted  that  zinc  has  a  negative  Poisson  ratio  vi2=-0.06.  In  Table  7  the  results  obtained 
by  Scholar  [13]  with  a  different  numerical  formulation  for  general  anisotropic  bodies  are  given  in 
the  second  column.  Form  a  comparison  between  columns  1  and  3  and  2  and  3,  it  turns  out  that 
the  present  formulation  and  implementation  is  much  more  precise. 


7.  Conclusions 

The  implementation  of  an  analytic  solution  for  Green’s  displacements  in  general 
anisotropic  solids  is  presented.  Its  detailed  illustration  has  been  accompanied  with  excerpts  from 
the  authors’  own  Fortran  code,  in  order  for  the  implementation  to  be  available  and  readily  usable 
by  as  many  readers  as  possible.  Many  features  distinguish  the  present  implementation  from  the 
existing  numerical  formulations: 

1)  the  procedure  is  completely  analytic,  the  only  numerical  step  is  represented  by  the 
determination  of  the  roots  of  a  sixth-order  polynomial; 

2)  once  the  roots  of  this  polynomial  are  known,  the  entire  Green’s  tensor  is  immediately 
calculated; 

3)  the  procedure  is  very  robust,  since  no  problem  arose  even  with  highly  degenerate  materials 
such  as  transversely  isotropic  or  isotropic  materials; 

4)  the  implementation  is  very  efficient,  since  less  than  16  sec  were  necessary  to  run  10,000 
calculations  of  the  entire  Green’s  tensor  in  a  PC  featuring  a  266  Mhz  Pentium  II  processor  and 
64  MB  RAM; 

5)  an  extensive  numerical  validation  (one  example  of  which  was  included  in  the  paper)  has 
shown  its  high  accuracy. 
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A  numerical  algorithm  has  also  been  proposed  for  the  Green’s  stresses  and  their 
derivatives.  Despite  its  simplicity,  it  has  been  proven  to  be; 

1)  robust,  since  no  problem  arose  even  with  highly  degenerate  materials  such  as  transversely 
isotropic  or  isotropic  materials; 

2)  very  accurate  even  with  degenerate  materials  and/or  when  the  field  point  is  very  close  to  or 
very  far  from  the  source  point ; 

3)  very  efficient,  since  less  than  80  sec  were  necessary  to  run  10.000  calculations  of  the  complete 
Green’s  stresses  in  a  PC  featuring  a  266  MHz  Pentium  II  processor  and  64  MB  RAM. 

Finally,  the  performance  of  the  proposed  implementations  within  a  previously  developed 
3D  BEM  code  [UJ  turned  out  to  be  highly  acerbate  when  compared  to  both  exact  solutions  and 
transversely  isotropic  BEM  formulations  for  which  the  closed-form  expressions  of  the  Green’s 
displacements,  stresses  and  stress  derivatives  were  used.  When  compared  to  previously  published 
results  obtained  with  completely  numerical  formulations,  the  present  implementation  turned  out 
to  be  much  more  precise. 

In  conclusion,  the  anisotropic  Green’s  functions  and  the  corresponding  3D  BEM  code  are 
accurate  and  efficient.  The  authors  are  currently  investigating  the  effect  of  rock  anisotropy, 
irregular  topography,  and  underground  openings  on  the  3D  in-situ  stresses.  The  outcome  of  such 
an  analysis  could  have  application  in  some  of  the  air  force  projects,  in  particular,  in  the 
penetration  mechanics  related  project. 
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10.  Appendix  A:  Radon  Transform  and  plane  representation  of  the  Dirac 
delta  function 

The  Radon  Transform  ([14]-[19])  is  of  fundamental  importance  in  order  to  work  out  the 
analytic  solution  of  the  problem  at  hand. 

Let  Ax)  be  a  function  defined  in  and  s  a  real  number;  the  Radon  transform  of  Ax)  is 
defined  as : 

f{s,n)  =  R|/(x)]  =  J/(x)-  <^(5  -n  •  x)dx  (A.  1) 

where  80  is  the  one-dimensional  Dirac  delta  function. 

It  follows  that,  when  s  varies  over  the  real  line,  the  Radon  transform  is  an  integration  of 
Ax)  over  all  planes  defined  by  n  •  x  =  j ,  i.e.  having  normal  n  and  distant  s\  n  I  from  the  origin  O. 

The  inverse  Radon  transform  is  an  integration  in  the  n-space  over  the  closed  surface  Q 
containing  the  origin,  defined  as  : 

f(x)=R  •  (/■■)=  f/"(n ■  i,ii>in(n)  (A.2) 

<S7I  ' 


where : 


ds^ 


(A.3) 


Let  <^(x)=  S{x^,X2,x^  be  the  Dirac  delta  centered  in  the  origin,  i.e.  the  functional : 
J^(x)/(xVF  =  /(o)  (A.4) 

where  o=(0,0,0). 

We  will  use  the  same  symbol  6  for  both  one-dimensional  and  three-dimensional  Dirac  delta,  with 
the  convention  that  if  the  argument  is  a  scalar,  the  one-dimensional  Dirac  delta  is  involved  and  if 
the  argument  is  a  vector,  the  three-dimensional  Dirac  delta  is  involved. 

The  Radon  transform  of  the  Dirac  delta  is  : 


5{s,  n)  =  if  [<?(x)]  =  (x)'  -  n  •  x)d%.  =  ^(5  -  n  •  0)  = 


Now: 


5<5(n-x)_  5^<y(n  x)_  2d^S 

ds  ’  dx,^  d^s 


(A.5) 


(A.6.a,b) 


Thus : 

•^5^^(n-x)  2d^S 

,=i  dXf  ds 

Since  the  first  member  of  (A.7)  is  Laplacian  of  S,  AS,  Eq.  (A.3)  becomes: 

d^d\  _A5(n-x) 


_df5 

2  1  |2  5 

ds^ 

5=n-x 

5=n'X 

(A.7) 


d"=- 


ds^ 


(A.8) 
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According  to  (A.2)  the  inverse  Radon  transform  is  : 

>  f^£%4«(n)=  -^A 

'  ’  J  l|2  ^  8;r  •*  ^  ’ 


n  n 


(A.9) 


n  n 


the  last  passage  is  due  to  the  fact  that  the  variable  of  integration  is  n,  not  x. 

Thus,  we  have  the  very  notable  relation,  called  “plane  representation  for  ^x)”: 

j(x)  =  -^Aj^^j2^Q(n)  (A.  10) 

8^  o  |ii| 


that  coincides  with  Eq.  (6)  in  Wang  [3]. 


11.  Appendix  B:  Change  of  coordinate  system 

The  Radon  transform  is  an  integration  over  the  planes  whose  normal  is  n.  The  inverse 
Radon  transform,  for  a  fixed  x,  is  an  integration  involving  all  the  normal  vectors  n.  Therefore,  a 
convenient  coordinate  system  when  we  perform  the  inverse  transform  is  such  that  an  axis  is 
parallel  to  x  [3]  (see  Figure  l.a).  Let  us  define: 

r  =  |x|;  e  =  ^  (B.l) 

If  v  is  an  arbitrary  unit  vector  different  fi-om  e  (v^^e),  two  normal  vectors  orthogonal  to  e 

are : 

ex  v 

P  =  l - r  (B.2) 

q  =  exp  (B.3) 

Let  4,  (^,  q  be  the  components  of  vector  n  in  the  new  coordinate  system  of  then: 
n  =  4P  +  ^q  +  7e  (B.4a) 

n-x  =  p-x^  +  q-xC  +  e-x;7  =  r7  (B.4b) 

This  transformation  induces  a  transformation  of  coordinates  in  the  n-space  with  (^,  (;,  r|) 
being  the  coordinates  of  point  n  in  the  new  coordinate  system  in  n-space,  see  Figure  1  .b.  The 
determinant  of  the  Jacobian  of  the  latter  transformation  is  obviously  equal  to  1 . 


12.  Appendix  C:  coefficients  byk 

If  v=(l,0,0)  then: 

b(l,  1,  1)  =  (yf**2*cd(5,5)  -  2*yf*zf*cci(5, 6)  +  zf**2*cd (6, 6) )  / 
-(yf**2  +  zf**2) 

C 

b(l,  1,  2)  =  2* (yf**3*cd(l,5)  +  yf*zf**2*cd(l, 5) -yf**2*zf*cd(l, 6) - 
-zf**3*cd(l, 6)  -  xf*yf*zf*cd(5,5)  ~  xf *yf**2*cd (5, 6)  + 
-xf*zf**2*cd(5,6)  +  xf*yf*zf*cd(6,6) )/ 

-((yf**2  +  zf**2) *Sqrt {xf**2  +  yf**2  +  zf**2) ) 

C 

b(l,  1,  3)  =  (yf**4*cd(l,l)  +2*yf**2*zf**2*cd(l,l)+zf**4*cd(l,l)- 
-2*xf*yf**2*zf*cd(l,5)-2*xf*zf**3*cd(l,5)-2*xf*yf**3*cd(l, 6)- 
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-2*xf*yf*zf**2*cd(l, 6)  +  xf**2*zf**2*cd(5, 5) +2*xf**2*yf*zf*cd(5, 6)+ 
-xf**2*yf**2*cd(6,6))/{(yf**2  +  zf**2) * (xf**2  +  yf**2  +  zf**2)) 

C 

b{l,  2,  1)  -  (“(yf*zf*cd(2,5) )  +  zf **2*cd (2, 6)  +  yf **2*cd (4, 5)  ^ 
-yf*zf*cd(4,6) )/(yf**2  +  zf**2) 

C 

b(l,  2,  2)  -  {-(yf**2*zf*cd{l,2)  )  -  zf**3*cd {1, 2 )  +  yf**3’"cd(l,  4)  + 
-yf*zf**2*cd(l,4)  -  xf*yf**2*cd(2,5)  +  xf*zf**2*cd(2, 5)  + 

~2*xf *yf*zf *cd (2, 6)  -  2*xf*yf*zf *cd {4, 5)  -  xf*yf **2*cd (4, 6)  + 
-xf*zf**2*cd(4,6)  +  yf**3*cd(5,6)  +  yf*zf**2*cd(5,6)“ 
-yf**2*zf*cd(6,6)  -  zf **3*cd ( 6, 6) ) / 

-((yf**2  +  zf**2) ♦Sqrt (xf **2  +  yf**2  +  zf**2) ) 

C 

b(l,  2,  3)  =  (“(xf*yf**3*cd(l,2) )  -  xf*yf*zf **2*cd {1, 2)  - 
-xf*yf**2*zf *cd (1, 4)  -  xf*zf**3*cd(l, 4)  +  yf**4*cd{l, 6)  + 
-2*yf**2*zf**2*cd(l,6)  +  zf**4*cd (1, 6)  +  xf**2*yf *zf*cd (2, 5)  + 
''Xf**2*yf**2^t5d{2, 6)  +  xf **2*zf**2^cd (4, 5)  ^  xf**2*yf *zf *cd (4,  6)  - 
-xf*yf**2*zf*cd(5, 6)  -  xf*zf**3*cd(5,6)  -  xf *yf **3*cd (6, 6)  - 
-*xf*yf*zf**2*cd(6,6))/((yf**2  +  zf**2)  *  (xf**2  +  yf**2  +  zf**2)) 

C 

b(l,  3,  1)  -  (yf**2*cd(3,5)  -  yf *zf*cd (3, 6)  ”  yf *zf*cd (4, 5)  + 
-zf**2*cd(4,6))/(yf**2  +  zf**2) 

C 

b(l,  3,  2)  =  (yf**3*cd(l,3)  +  yf  *zf  **2*cd  (1, 3)  -  yf  **2*zf  *cd  (1,  4) -> 
-zf**3*cd(l,4)  -  2*xf*yf*zf*cd(3, 5)  -  xf*yf**2*cd(3, 6)  + 
-xf*zf**2*cd(3,6)  -  xf*yf**2*cd(4,5)  +  xf*zf**2*cd(4, 5)  + 
-2*xf*yf*zf*cd(4,6)  +  yf**3*cd(5, 5)  +  yf*2f**2*cd (5, 5)  - 
-yf**2*zf*cd(5,6)  -  zf **3*cd (5, 6) ) / 

-((yf**2  +  zf**2)*Sqrt{xf**2  +  yf**2  +  zf**2)) 

C 

b(l,  3,  3)  =  (-(xf*yf**2*zf*cd(l,3) )  -  xf*zf**3*cd(l,3)- 
-xf*yf**3*cd(l,4)  - 

-xf*yf*zf**2*cd(l,4)  +  yf**4*cd (1, 5)  +  2*yf**2*zf**2*cd (1, 5)  + 
-‘Zf**4*cd(l,5)  +  xf**2*zf**2*cd(3,5)  +  xf**2*yf *zf*cd (3, 6)  + 
-xf**2*yf*zf*cd(4,5)  +  xf **2*yf**2*cd(4, 6)  -  xf*yf**2*zf*cd (5, 5)  - 
-xf*zf**3*cd(5,5)  -  xf*yf**3*cd(5,6)  -  xf ^yf *zf **2*cd (5, 6) ) / 
~((yf**2  +  zf**2) * (xf**2  +  yf**2  +  zf**2)) 

C 

b(2,  1,  1)  -  b(l,  2,  1) 

b{2,  1,  2)  =  b(l,  2,  2) 

b(2,  1,  3)  -  b(l,  2,  3) 

C 

b(2,  2,  1)  =  (zf**2*cd(2,2)  -  2*yf *zf*cd {2, 4 )  +  yf **2*cd (4, 4) ) / 
-(yf**2  +  zf**2) 

C 

b(2,  2,  2)  =  2*(xf*yf*zf*cd(2,2)  -  xf*yf**2*cd (2, 4)  + 
-xf*zf**2*cd(2,4)  - 

-yf**2*zf*cd(2,6)  -  zf**3*cd(2,6)  ~xf*yf*2f*cd(4, 4)  +yf**3*cd(4,6)+ 
-yf*zf**2*cd(4,6))/((yf**2  +  zf**2) *Sqrt (xf **2  +  yf**2  +  zf**2) ) 

C 

b(2,  2,  3)  =  {xf**2*yf**2*cd(2,2)  +  2*xf**2*yf*2f*cd(2, 4)  ~ 
-2*xf*yf**3*cd{2,6)  -  2*xf *yf *zf**2*cd(2, 6)  +  xf**2*zf**2*cd (4, 4)  - 
-2*xf*yf**2*2f*cd{4,6)  -  2*xf*zf **3*cd (4, 6)  +  yf**4*cd(6, 6)  + 
“2*yf**2*zf**2*cd(6,6)  +  zf**4*cd{6, 6) ) / 

-’{(yf**2  +  zf**2)*(xf**2  +  yf**2  +  zf**2)) 

C 

b(2,  3,  1)  -  (~(yf*zf*cd(2,3)  )  +  zf**2*^cd(2, 4)  +  yf**2*cd(3, 4)  - 
-yf*zf*cd(4,4) )/(yf**2  +  zf**2) 

C 

b(2,  3,  2)  =  (-(xf*yf**2*cd(2,3) )  +  xf*zf**2*cd(2,3)  + 
-2*xf*yf*zf*cd(2,4)  - 

-yf**2*zf*cd(2,5)  -  zf**3*cd(2,5)-2*xf*yf*zf*cd(3,4)+yf**3*cd(3, 6)  + 
-yf*zf**2*cd(3,6)  -  xf*yf**2*cd(4, 4)  +  xf*zf**2*cd (4, 4 )  + 
“yf**3*cd(4,5)  +  yf *2f**2*cd (4, 5)  -yf**2*zf*cd (4, 6) “2f**3*cd (4, 6) ) / 
-((yf**2  +  zf**2)*Sqrt{xf**2  +  yf**2  +  zf**2)) 

C 

b(2,  3,  3)  =  (xf**2*yf*zf*cd(2,3)  +  xf**2*yf**2*cd (2, 4 )  - 
“Xf*yf**3*cd(2,5)  -* 

-xf*yf*zf**2*cd(2,5)  +  xf**2*zf**2*cd (3, 4)  -  xf*yf **2*zf*cd (3, 6)  - 
-xf*zf**3*cd(3, 6)  +  xf**2*yf*zf*cd(4,4)  -  xf*yf**2*zf *cd (4, 5)  - 
-xf*zf**3*cd(4,5)  -  xf*yf**3*cd(4, 6)  -  xf*yf*zf**2*cd (4, 6)  + 
-yf**4*cd(5,6)  +  2*yf**2*zf**2*cd(5,6)  +  zf**4*cd (5, 6) ) / 

-((yf**2  +  zf**2) * (xf**2  +  yf**2  +  zf**2) ) 


b(3. 

1,  1)  -  b(l, 

3, 

1) 

b(3. 

1,  2)  =  b(l. 

3, 

2) 

b(3. 

1,  3)  =  b(l. 

3, 

3) 

b(3. 

2,  1)  =  b(2, 

3, 

1) 

b(3. 

2,  2)  »  b(2, 

3, 

2) 
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c 


c 


c 


b(3,  2,  3)  *  b(2,  3,  3) 

b{3,  3,  1)  =  (yf**2*cd(3,3)  -  2*yf*zf *cd (3, 4 )  +  zf **2*cd {4, 4 ) } / 
“(yf**2  +  zf**2) 

b(3,  3,  2)  -  2*{-(xf*yf*zf*cd(3,3) )  -  xf *yf **2*cd (3, 4 )  + 
-xf*zf**2*cd(3,4)  + 

-yf**3*cd (3, 5)  +  yf*zf **2*cd {3, 5)  +  xf*yf*zf*cd(4, 4)  - 
-yf**2*zf*cd(4,5)  -  zf**3*cd(4, 5) ) / 

-((yf**2  +  zf*^2)*Sqrt(xf**2  +  yf*'^2  +  zf^*2)  ) 

b(3,  3,  3)  -  (xf**2*zf**2*cd(3,3)  +  2*xf **2*yf*zf*cd (3, 4)  - 
-2*xf*yf**2*zf*cd(3,5)  -  2*xf *zf**3*cd (3, 5)  +  xf **2*yf**2*cd (4, 4) 
~2*xf*yf**3*cd(4, 5)  -  2*xf*yf*zf**2*cd (4, 5)  +  yf **4*cd (5, 5)  + 
-2*yf**2*zf**2*cd(5,5)  +  zf **4*cd (5, 5) ) / 

-((yf**2  +  2f**2)*(xf**2  +  yf**2  +  zf**2) ) 


13,  Appendix  D:  coefficients  ai 

a(l)=-(b{l,3,l)*b(2,2,l)*b(3,l,l))  + 
-2*b(l,2,l)*b(2,3,l)*b(3,l,l)- 
-b(l,l,l)*b(2,3,l)*b(3,2,l)- 
-b(l,2,l)*b(2,l,l)*b(3,3,l)  + 
-b(l,l,l)*b(2,2,l)*b(3,3,l); 
c 

a(2)=-(b(l,3,2)*b(2,2,l)*b(3,l,l)  + 
-b(l,3,l)*b(2,2,2)*b(3,l,l)  + 
-b{l,3,l)*b(2,2,l)*b(3,l,2))  + 
-2*(b(l,2,2)*b(2,3,l)*b(3,l,l)  + 
-b(l,2,l)*b(2,3,2)*b(3,l,l)  + 
-b(l,2,l)*b(2,3,l)»b{3,l,2))- 
-(b(l,l,2)*b(2,3,l)*b(3,2,l)  + 
-b(l,l,l)*b(2,3,2)*b(3,2,l)  + 
-b(l,l,l)*b(2,3,l)*b(3,2,2))- 
-(b(l,2,2)*b(2,l,l)*b(3,3,l)  + 
-b(l,2,l)*b(2,l,2)*b(3,3,l)  + 

-b(l,2,l)*b(2,l,l)*b(3,3,2))  + 
-(b(l,l,2)*b(2,2,l)*b(3,3,l)  + 
-b(l,l,l)*b(2,2,2)*b(3,3,l)  + 

-b(l,  1,1)  *b(2, 2,1)  *b (3,3,2)  ); 

c 

a(3)=-(b(l,3,3)*b(2,2,l)*b(3,l,l)  + 
-b(l,3,l)*b(2,2,3)*b(3,l,l)  + 
-b(l,3,l)*b(2,2,l)*b(3,l,3)  + 
-b(l,3,2)*b(2,2,2)*b(3,l,l)  + 
-b(l,3,l)*b(2,2,2)*b(3,l,2)  + 
-b(l,3,2)*b(2,2,l)*b(3,l,2))  + 
-2*(b(l,2,3)*b(2,3,l)*b(3,l,l)  + 
-b(l,2,l)*b(2,3,3)*b(3,l,l)  + 
-b(l,2,l)*b(2,3,l)*b(3,l,3)  + 
-b(l,2,2)*b(2,3,2)*b(3,l,l)  + 
-b(l,2,l)*b(2,3,2)*b(3,l,2)  + 
-b(l,2,2)*b(2,3,l)*b(3,l,2))- 
-(b(l,l,3)*b(2,3,l)*b(3,2,l)  + 
-b(l,l,l)*b(2,3,3)*b(3,2,l)  + 
-b(l,l,l)*b(2,3,l)*b(3,2,3)  + 
-b(l,l,2)*b(2,3,2)*b(3,2,l)  + 
-b(l,l,l)*b(2,3,2)*b(3,2,2)  + 
-b(l,l,2)*b(2,3,l)*b(3,2,2))- 
-(b(l,2,3)*b(2,l,l)*b(3,3,l)  + 
-b(l,2,l)*b(2,l,3)*b(3,3,l)  + 
-b(l,2,l)*b(2,l,l)*b(3,3,3)  + 
-b(l,2,2)*b(2,l,2)*b(3,3,l)  + 
-b(l,2,l)*b(2,l,2)*b(3,3,2)  + 
-b(l,2,2)*b(2,l,l)*b(3,3,2)  )  + 
-(b(l,l,3)*b(2,2,l)*b{3,3,l)  + 
-b(l,l,l)*b(2,2,3)*b(3,3,l)  + 
-b(l,l,l)*b(2,2,l)*b(3,3,3)  + 
-b(l,l,2)*b(2,2,2)*b(3,3,l)  + 
-b(l,l,l)*b(2,2,2)*b(3,3,2)  + 
-b(l,l,2)*b(2,2,l)*b(3,3,2)) 

a(4)=-(b(l,3,3)*b(2,2,2)*b(3,l,l)  + 
-b(l,3,2)*b(2,2,3)*b(3,l,l)  + 
-b(l,3,l)*b(2,2,2)*b(3,l,3)  + 
-b(l,3,l)*b(2,2,3)*b(3,l,2)  + 
-b(l,3,2)*b(2,2,l)*b(3,l,3)  + 
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-b{l,3,3)*b{2,2,l)*b(3,l,2)  + 
-b(l,3,2)*b(2,2,2)*b(3,l,2))  + 
-2*(b(l,2,3)*b(2,3,2)*b(3,l,l)  + 
-b(l,2,2)*b(2,3,3)*b(3,l,l)  + 
-b{l,2,l)*b(2,3,2)*b(3,l,3)  + 
-b{l,2,l)*b(2,3,3)*b(3,l,2)  + 
-'b(l,2,2)*b(2,3,l)*b(3,l,3)  + 
-b(l,2,3)*b(2,3,l)*b(3,l,2)  + 
-b(l,2,2)*b(2,3,2)*b(3,l,2))- 
-(b{l,l,3)*b(2,3,2)*b{3,2,l)  + 
-b(l,l,2)*b(2,3,3)*b(3,2,l)  + 
-b(l,l,l)*b(2,3,2)*b(3,2,3)  + 
-b(l,l,l)*b(2,3,3)*b(3,2,2)  + 
-b(l,l,2)*b(2,3,l)*b(3,2,3)  + 
-b(l,l,3)*b(2,3,l)*b(3,2,2)  + 
-b(l,l,2)*b(2,3,2)*b(3,2,2))- 
-(b(l,2,3)*b(,2^1,2)*b(3,3,l)  + 
-b(l,2,2)*b(2,l,3)*b{3,3,l)  + 
“b{l,2,l)*b{2,l,2)*b(3,3,3)+ 
-b(l,2,l)*b(2,l,3)*b(3,3,2)  + 
-b(l,2,2)*b(2,l,l)*b(3,3,3)  + 
-b(l,2,3)*b(2,l,l)*b(3,3,2)  + 
-b{l,2,2)*b(2,l,2)*b(3,3,2))+ 
-(b(l,l,3)*b(2,2,2)*b(3,3,l)  + 
-b(l,l,2)*b(2,2,3)*b(3,3,l)  + 
-b{l,l,l)*b(2,2,2)*b(3,3,3)  + 
-b(l,l,l)*b(2,2,3)*b(3,3,2)  + 
-b(l,l,2)*b(2,2,l)*b(3,3,3)  + 
-b(l,l,3)*b(2,2,l)*b(3,3,2)  + 
-b(l,l,2)*b(2,2,2)*b(3,3,2)) 

a(5)«-(b{l,3,2)*b{2,2,2)*b(3,l,3)  + 
-b(l,3,3)*b(2,2,2)*b(3,l,2)  + 
-b(l,3,2)*b(2,2,3)*b{3,l,2)  + 
-b(l,3,l)*b(2,2,3)*b(3,l,3)  + 
-b(l,3,3)*b(2,2,l)*b(3,l,3)  + 
-b(l,3,3)*b(2,2,3)*b(3,l,l))+ 
-2*(b(l,2,2)*b(2,3,2)*b(3,l,3)+ 
-b{l,2,3)*b(2,3,2)*b(3,l,2)  + 
-b{l,2,2)*b(2,3,3)*b(3,l,2)  + 
-b(l,2,l)*b(2,3,3)  *b(3, 1,3)  + 
-b(l,2,3)*b(2,3,l)*b(3,l,3)  + 
-b(l,2,3)*b(2,3,3)*b(3,l,l))- 
-(b(l,l,2)*b(2,3,2)*b(3,2,3)  + 
-b{l,l,3)*b(2,3,2)*b(3,2,2)  + 
-b{l,l,2)*b(2,3,3)*b(3,2,2)  + 
-b(l,l,l)*b(2,3,3)*b(3,2,3)  + 
“b(l,l,3)*b(2,3,l)*b(3,2,3)  + 
-b(l,l,3)*b(2,3,3)*b{3,2,l)>- 
-(b(l,2,2)*b(2,l,2)*b(3,3,3)  + 
-b{l,2,3)*b(2,l,2)*b(3,3,2)  + 
-b(l,2,2)*b(2,l,3)*b(3,3,2)  + 
-b(l,2,l)*b(2,l,3)*b(3,3,3)  + 
-b(l,2,3)*b(2,l,l)*b(3,3,3)  + 
-b(l,2,3)*b(2,l,3)*b(3,3,l))+ 
-{b(l,l,2)*b(2,2,2)*b(3,3,3)+ 
-b(l,l,3)*b(2,2,2)*b(3,3,2)  + 
-b(l,l,2)*b(2,2,3)*b(3,3,2)  + 
-b(l,l,l)*b(2,2,3)*b(3,3,3)  + 
-b(l,l,3)*b(2,2,l)*b(3,3,3)  + 
-b(l,l,3)*b(2,2,3)*b(3,3,l)  ) 

a{6)=-(b(l,3,2)*b(2,2,3)*b(3,l,3)  + 
~b(l,3,3)*b(2,2,2)*b(3,l,3)  + 
-b(l,3,3)*b(2,2,3)*b(3,l,2))  + 
-2*(b{l,2,2)*b(2,3,3)*b(3,l,3)  + 
-b|l,2,3)*b(2,3,2)*b(3,l,3)  + 
-b(l,2,3)*b(2,3,3)*b(3,l,2))- 
-(b(l,l,2)*b(2,3,3)*b(3,2,3)+ 
-b(l,l,3)*b(2,3,2)*b(3,2,3)  + 
-b(l,l,3)*b(2,3,3)*b(3,2,2))- 
-(b(l,2,2)*b(2,l,3)*b(3,3,3)  + 
-b(l,2,3)*b(2,l,2)*b(3,3,3)  + 
-b(l,2,3)*b(2,l,3)*b(3,3,2))  + 
-(b(l,l,2)*b(2,2,3)*b(3,3,3)  + 
-b(l,l,3)*b(2,2,2)*b(3,3,3)+ 
-b(l,l,3)*b(2,2,3)*b(3,3,2)  ) 


a(7)=-(b(l,3,3)*b(2,2,3)*b(3,l,3))  + 


-2*b(l,2,3)*b(2,3,3)*b(3,l,3)- 
-b(l,l/3)*b(2,3,3)*b(3,2,3)- 
-b(l,2,3)*b(2,l,3)*b(3,3,3)  + 
-b(l,l,3)*b{2,2,3)*b{3,3,3) 


Figure  captions: 

Fig.  l.a.  Anisotropic  elastic  body  (shown  bounded  for  representation  convenience)  referred  to  a 
fixed  Cartesian  system. 

Fig.  l.b.  Closed  surface  Q  containing  the  origin  in  the  n-space. 

Fig.  2.  Parallelepiped  over  which  contour  integration  is  taken  in  the  n-space;  Surfaces  Si  and  S2 
are  boimded  by  ^=±1,  S3  and  S4  by  ^=±L,  S5  and  Ss  by  ri=±L. 

Fig.  3.  Cube  having  edge  of  length  /=!  m  discretized  with  six  nine-node  quadrilateral  elements 
with  a  total  of  26  nodes. 

Fig.  4.  Orientation  of  the  plane  of  transverse  isotropy.  Dip  angle  \}/  is  the  angle  between  the  plane 
of  symmetry  and  xiX2  plane;  dip  direction  angle  9  is  the  angle  between  X2  and  the  orthogonal 
projection  of  the  dip  vector  on  the  X1X2  plane.  Positive  angles  are  shown. 
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Fig.  3.  Cube  having  edge  of  length  /=1  m  discretized  with  six  nine-node  quadrilateral  elements 
with  a  total  of  26  nodes.  Nodes  from  1  to  8  have  fixed  displacement  along  X3,  node  9  has  fixed 

displacements  along  jci  and  X2. 
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PLANE  OF 

TRANSVERSE 

ISOTROPY 


DIP  DIRECTION 


DIP  VECTOR 


Fig.  4.  Orientation  of  the  plane  of  transverse  isotropy.  Dip  angle  is  the  angle  between  the  plane 
of  symmetry  and  xiX2  plane;  dip  direction  angle  (p  is  the  angle  between  X2  and  the  orthogonal 
projection  of  the  dip  vector  on  the  xixz  plane.  Positive  angles  are  shown. 
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Table  1.  Correspondence  between  notation  used  in  the  text  and  notation  used  in  the  present 


FORTRAN  code. 


Notation  used  in  the  text 

Notation  used  in  the  Fortran  code 

biik 

b(ij,k) 

_ Dii 

cd(ij) 

X=(Xi,  X2,  X3) 

xf,  yf,  zf 

ai  a(i) 

Table  2.  Greenes  displacements  (m"^)  calculated  according  to  Pan’s  and  Chou’s  closed-form 
solution  [4]  as  implemented  in  [1 1]  and  with  the  present  formulation.  The  source  point  is  at  the 


origin,  the  field  point  is  at  x=(- 1,0.8, 1.5)  m. 


as 

- — D— ?  — - r-- 

gij  transversely 
isotropic 
formulation 

gij  present 
formulation 

relative  difference 

1,1 

4.0141588565E-03 

4.0141588610E-03 

l.lE-09 

1,2 

-2.9315529284E-04 

-2,9315529143E-04 

4.8E-09 

1,3 

-2.1517885087E-03 

-2.1517885172E-03 

3. 9E-09 

2,1 

-2.9315529284E-04 

-2.9315529143E-04 

4 .8E-09 

2,2 

3.8822389747E-03 

3.8822389799E-03 

1.3E-09 

2,3 

1.7214308070E-03 

1.7214308137E-03 

3.9E-09 

3,1 

-2.1517885087E-03 

-2'.1517885172E-03 

3.9E-09 

3,2 

1.7214308070E-03 

1.7214308137E-03 

3.9E-09 

3,3 

1. 9003220124E-02 

1.9003220284E-02 

8.4E-09 
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Table  3.  Green’s  stresses  (kN/m^)  calculated  according  to  Pan’s  and  Chou’s  closed-form 
solution  [4]  as  implemented  in  [1 1]  and  the  present  formulation.  The  source  point  is  at  the  origin, 
_ the  field  point  in  x=(-l»0  8,1.5)  m. _ 


ijk 

Oijk  transversely 
isotropic 
formulation 

Gijk  present 
formulation 

relative 

difference 

111 

3.0534645334E-03 

3.0534575287E-03 

2.3E-06 

221 

5.5516532436E-03 

5.5516460542E-03 

1.3E-06 

331 

6.9594542712E-03 

6,9594505349E-03 

5.4E-07 

231^^^ 

2.5160305656E-0fi 

2.5160307422E-03 

7  .OE-08 

131 

-3.3933330101E~03 

-3,3933332396E-03 

6.8E-08 

121 

-1.8033628216E-03 

-1.8033626725E-03 

8.3E~08 

111 

-5.6717491681E-03 

-5.6717487868E-03 

6.7E~08 

221 

"1.2123450534E-03 

~1.2123450967E-03 

3.6E-'08 

331 

-5.5675634169E-03 

-5.5675634915E-03 

1.3E-08 

231 

-2,2611192556E-03 

-2.2611192104E-03 

2.0E-08 

131 

2.5160305656E-03 

2.5160306251E-03 

2.4E~08 

121 

7.1617031049E-04 

7.1617003567E-04 

3.8E-07 

111 

-1.7147827999E-02 

-1.7147807646E-02 

1.2E-06 

221 

-8.3655438558E-03 

-8.3655210747E-03 

2.7E-06 

331 

-2.2078536037E-02 

-2.2078524002E-02 

5.4E-07 

231 

-1.1775219220E-02 

-1.1775219467E-02 

2,lE-08 

131 

1.4719024025E-02 

1.4719024871E-02 

5.7E-08 

121 

1.9516186984E-02 

1.9516189072E-02 

l.lE-07 
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Table  4.  Derivatives  of  the  Green’s  stresses  (kN/m^)  calculated  according  to  Pan’s  and  Chou’s 
closed-form  solution  [4]  as  implemented  in  [1 1]  and  the  present  formulation.  The  source  point  is 


at  the  origin,  the  field  point  in  x=(-l»0-8,l  .5)  m. 


Ijk,l 

tJijk,!  transversely 
isotropic 
formulation 

<7ijk,i  present 
formulation 

relative 

difference 

111,1 

-0,59622259271E-02 

-0.59629572705E-02 

1.2E-04 

111,2 

~0.45086940696E“02 

>-0.45074656599E-02 

“2.7E-04 

111,3 

~0.56414664921E-02 

-'0.56411267186E-02 

-6.0E-05 

221,1 

0.44694421882E-02 

0.44685063303E-02 

-2.0E-04 

221,2 

-0.11811731607E“02 

-0.11798500052E-02 

-l.lE-03 

221,3 

‘-0.37926171802E-fi2 

. . .  ii  iiimimm 

“1.3E-04 

331,1 

-3.2E-04 

331,2 

-0.69265941556E-02 

-0.69258768984E-02 

-l.OE-04 

331,3 

-0.44525631963E-02 

-0.44523121121E-02 

-5.6E-05 

231,1 

0.79385810234E-03 

0.79379727638E-03 

-7.7E-05 

231,2 

0.49712727263E-03 

0.49714271358E-03 

3.1E-05 

231,3 

-0.30906032310E-02 

-0.30906216080E-02 

5.9E-06 

131,1 

0.39554359236E-02 

0.39554842055E-02 

1.2E-05 

131,2 

0.18677123922E-02 

0.18676962357E-02 

“8.6E-06 

131,3 

0.61652880201E-02 

0,61653303015E-02 

6.8E-06 

121,1 

0.42717763920E-02 

0.42718112393E-02 

8.1E-06 

121,2 

-0.20306209291E-03 

-0.20300240322E-03 

-2.9E-04 

121,3 

0.53606344730E-02 

0.53607252836E-02 

1.7E-05 

112,1 

-0.30768883497E-02 

-0.30775369921E-02 

2.1E-04 

112,2 

0.84038676742E-03 

0.84094663232E-03 

6.7E-04 

112,3 

0.50628670483E-02 

0.50629208004E-02 

1.6E-05 

222,1 

-0.26129788805E-02 

-0.26137365442E-02 

2.9E-04 

222,2 

-0.48936107600E-02 

-0.48931499059E-02 

-9.4E-05 

222,3 

0.24843998895E-02 

0.24845314971E-02 

5.3E-05 

332,1 

~0.69265941556E-02 

-0.69269468094E-02 

5.9E-05 

332,2 

-0.14181789466E-02 

-0.14179740950E-02 

-l,4E-04 

332,3 

0.35620505570E-02 

0.35621376993E-02 

2.4E-05 

232,1 

-0.84519058971E-03 

-’0.84518000388E-03 

-1.2E-05 

232,2 

2.9E-06 

232,3 

0.47745165660E-02 

0.47746037165E-02 

1.8E-05 

132,1 

0.79385810234E~03 

0.79378559587E-03 

~9.lE-05 

132,2 

0.49712727263E-03 

0.49718754586E-03 

1.2E-04 

132,3 

-0.30906032312E-02 

-0.30906169537E-02 

4.4E-06 

122,1 

0.11909419404E-03 

0.11905631360E-03 

-3.2E-04 

122,2 

0.61674915810E-02 

0.61674862040E-02 

-8.7E-07 

122,3 

-0.41648264611E-02 

-0.41648490465E-02 

5.4E-06 

113,1 

0.71134008460E-02 

0.71131995957E-02 

-2.8E-05 

113,2 

3.4E-05 

113,3 

0.98238228068E-02 

0.98259170804E-02 

2.1E-04 

223,1 

-0.29318534460E-01 

-0.29319149219E-01 

3.0E-05 

223,2 

“0.15577546399E-01 

-0.15576273564E-01 

-8.2E-05 

223,3 

“0.83606419668E-04 

-0.81186914690E-04 

-2.9E-02 

333,1 

“0.27302017639E-01 

-0.27302079852E-01 

2.3E-06 

333,2 

0.21841614111E-01 

0.21842236628E-01 

2 . 8E-05 

333,3 

”0.41215790310E~03 

-0.41093285237E-03 

-3.0E-03 

233,1 

-0,14561076074E-01 

-0.14560962479E-01 

-7.8E~06 

233,2 

-0.30701631651E-02 

~0.30703156075E-02 

5.0E-05 

233,3 

0.76303285982E-02 

0.76304202029E-'02 

1.2E-05 

133,1 

0.34823210683E-02 

0.34819436621E-02 

-l.lE-04 

133,2 

-0.14561076074E-01 

-0.14560935326E-01 

-9.7E-06 

131,3 

-0.95379107477E-02 

“0.95379106793E-02 

-7.2E-09 

123,1 

0.79472178014E-02 

0.79465584372E-02 

-8.3E-05 

123,2 

0.24245099017E-02 

0.24250012946E-02 

2.3E-04 

123,3 

-0.22016509392E-01 

-0.22016463462E-01 

-2.1E-06 
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Table  5.  Block  under  uniform  compression.  Example  1 : 

_ C _ J* _ 1 _ A _ X _ 1-^  A 


Present 

Transv.  isotropic 
form,  fill 

Exact 

[12] 

Wi 

-0.5496 

-0.5496 

-0.5496 

U2 

-0.7029 

-0.7032 

-0.7031 

W3 

-1.0417 

-1.0417 

-1.0417 

Table  6.1.  Block  under  uniform  compression.  Example  1:  internal  displacements  along  the 

vertical  center  line. 


z(m) 

Present 

M3  (xlO’^  m) 

Transv.  isotropic  form.  [11] 

Exact  [12] 

0.75 

-0.5210 

-0.5208 

-0.5209 

0.625 

-0.2604 

-0.2604 

-0.2604 

0.5 

0.0000 

0.0000 

0.0000 

0.375 

0.2604 

0.2604 

0.2604 

0.25 

0.5208 

0.5208 

0.5208 

Table  6.2.  Block  xmder  uniform  compression.  Example  1: 
_ stresses  along  the  vertical  center  line. _ 


z(m) 

Present 

cii  (kN/m^) 

Transv.  isotropic  form.  [11] 

Exact  [12] 

0.75 

0.0005 

0.0011 

0.0000 

0.625 

0.0006 

-0.0000 

0.0000 

0.5 

-0.0023 

-0.0000 

0.0000 

0.375 

0.0096 

-0.0000 

0.0000 

0.25 

0.0013 

0.0014  . 

0.0000 
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Table  6.3.  Block  under  uniform  compression.  Example  1 ; 
_ stresses  along  the  vertical  center  line. _ 


z(m) 

Present 

C22  (kN/m^) 

Transv.  isotropic  form.  [Ill 

Exact  [121 

0.75 

-0.0060 

-0.0040 

0.0000 

0.625 

-0.0021 

-0.0000 

0.0000 

0.5 

-0.0007 

0.0000 

0.0000 

0.375 

0.0084 

-0.0000 

0.0000 

r.25 

-0.0045 

*-0.0044 

0.0000 

Table  6.4.  Block  under  imiform  compression.  Example  1: 
_ stresses  along  the  vertical  center  line. _ 


z(m) 

Present 

CT33  (kN/m^) 

Transv.  isotropic 
form.  [Ill 

Exact  [12] 

0.75 

-0.999 

-1.001 

-1.000 

0.625 

-0.992 

-1.000 

-1.000 

0.5 

-1.002 

-1.000 

-1.000 

0.375 

-1.000 

-1.000 

-1.000 

0.25 

-1.003 

0.988 

-1.000 

Table  7.  Zinc  block  under  uniform  compression.  Example  2: 

Kmin/tarv  at  no/lo  A  /"vl 


Present 

Anisotropic 
formulation  [131 

Exact 

[12] 

7.278 

7.291 

7.274 

U2 

-7.277 

-7.291 

-7.274 

W3 

-28.35 

-28.42 

-28.34 
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